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& ■ Abstract 

CN ■ Following a polynomial approach, many robust fixed-order controller design 

problems can be formulated as optimization problems whose set of feasible solutions 
Q^) | is modelled by parametrized polynomial matrix inequalities (PMI). These feasibility 

' sets are typically nonconvex. Given a parametrized PMI set, we provide a hierar- 

chy of linear matrix inequality (LMI) problems whose optimal solutions generate 
inner approximations modelled by a single polynomial superlevel set. Those inner 
approximations converge in a well-defined analytic sense to the nonconvex original 
feasible set, with asymptotically vanishing conservatism. One may also impose the 
hierarchy of inner approximations to be nested or convex. In the latter case they 
do not converge any more to the feasible set, but they can be used in a convex 
■ optimization framework at the price of some conservatism. Finally, we show that 

the specific geometry of nonconvex polynomial stability regions can be exploited to 
■^j- \ improve convergence of the hierarchy of inner approximations. 

Keywords: polynomial matrix inequality, linear matrix inequality, robust optimization, 
robust fixed-order controller design, moments, positive polynomials. 
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1 Introduction 

Linear system stability can be formulated semialgebraically in the space of coefficients 
of the characteristic polynomial. The region of stability is generally nonconvex in this 
space, and this is a major obstacle when solving fixed-order and/or robust controller 
design problems. Using the Hermite stability criterion, these problems can be formulated 
as parametrized polynomial matrix inequalities (PMIs) where parameters account for 
uncertainties and the decision variables are controller coefficients. Recent results on real 
algebraic geometry and generalized problems of moments can be used to build up a 
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hierarchy of convex linear matrix inequality (LMI) outer approximations of the region 
of stability, with asymptotic convergence to its convex hull, see e.g. [5] for a software 
implementation and examples, and see [7] for an application to PMI problems arising 
from static output feedback design. 

If outer approximations of nonconvex semialgebraic sets can be readily constructed with 
these LMI relaxations, inner approximations are much harder to obtain. However, for 
controller design purposes, inner approximations are essential since they correspond to 
sufficient conditions and hence guarantees of stability or robust stability. In the robust 
systems control literature, convex inner approximations of the stability region have been 
proposed in the form of polytopes [H], ellipsoids [4J or more general LMI regions [SJ [H] 
derived from polynomial positivity conditions. Interval analysis can also be used in this 
context, see e.g. [IS] . 

In this paper we provide a numerical scheme for approximating from inside the feasible 
set P C M" of a parametrized PMI P(x,u) y (for some matrix polynomial P), that is, 
the set of points x such that P(x, u) >z for all values of the parameter u in some specified 
domain Ucl p (assumed to be a basic compact semialgebraic set0). This includes as a 
special case the approximation of the stability region (and the robust stability region) of 
linear systems. The particular case where P(x, u) is affine in x covers parametrized LMIs 
with many applications in robust control, as surveyed e.g. in [T5] . 

Given a compact set B C W l containing P, this numerical scheme consists of building 
up a sequence of inner approximations G^ C P C B, d G N, which fulfils two essential 
conditions: 

1. The approximation converges in a well- defined analytic sense ; 

2. Each set G^ is defined in a simple manner, as a superlevel set of a single polynomial. 
In our mind, this feature is essential for a successful implementation in practical 
applications. 

More precisely, we provide a hierarchy of inner approximations (G^) of P, where each 
Gd = {x G B : gd(x) > 0} is a basic semi-algebraic set for some polynomial ga of degree d. 
The vector of coefficients of the polynomial is an optimal solution of an LMI problem. 
When d increases, the convergence of (G^) to P is very strong. Indeed, the Lebesgue 
volume of G^ converges to the Lebesgue volume of P. In fact, on any (a priori fixed) 
compact set B, the sequence (gd) converges for the Li-norm on B to the function x i— > 
Amin(^) = niin ue u A m i n (x, m) where A mm (x, w) is the minimum eigenvalue of the matrix- 
polynomial P(x, u) associated with the PMI. Consequently, gd —> A min in (Lebesgue) 
measure on B, and gd k — > A m i n almost everywhere and almost uniformly on B, for a 
subsequence (gd k )- In addition, if one defines the piecewise polynomial gd '■= max fc < rf g k , 
then cjd — > A m i n almost everywhere, almost uniformly and in (Lebesgue) measure on B. 

In addition, we can easily enforce that the inner approximations (G^) are nested and/or 
convex. Of course, for the latter convex approximations, convergence to P is lost if P is 

X A basic semialgebraic set is a set defined by intersecting a finite number of polynomial superlevel 
sets. 
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not convex. However, on the other hand, having a convex inner approximation of P may 
reveal to be very useful, e.g., for optimization purposes. 

On the practical and computational sides, the quality of the approximation of P depends 
heavily on the chosen set B D P on which to make the approximation of the function 
A m i n . The smaller B, the better the approximation. In particular, it is worth emphasizing 
that when the set P to approximate is the stability or robust stability region of a linear 
system, then its particular geometry can be exploited to construct a tight bounding set 
B. Therefore, a good approximation of P is obtained significantly faster than with an 
arbitrary set B containing P. 

Finally, let us insist that the main goal of the paper is to show that it is possible to 
provide a tight and explicit inner approximation with no quantifier, of nonconvex feasible 
sets described with quantifiers. Then this new feasible set can be used for optimization 
purposes and we are facing two cases: 

• the convex case: if / and — g are convex polynomials, B = {x 6 I™ : Halloo < 1} 
and G = {x G B : g{x) > 0} then the optimization problem min x f(x) s.t. x G 
G is polynomially solvable. Indeed, functions f(x), polynomially 
computable, of polynomial growth, and the feasible set is polynomially bounded. 
Then polynomial solvability of the problem follows from [31 Theorem 5.3.1]. 

• the nonconvex case: if — g is not convex then notice that firstly we still have an 
optimization problem with no quantifier, a nontrivial improvement. Secondly we 
are now faced with an polynomial optimization problem with a single polynomial 
constraint and possibly bound constraints x G B. One may then apply the hierarchy 
of convex LMI relaxations described in [11, Chapter 5]. Of course, in general, 
polynomial optimization is NP-hard. However, if the size of the problem is relatively 
small and the degree of g is small, practice seems to reveal that the problem is 
solved exactly with few relaxations in many cases, see pTT| §5.3.3]. In addition, if 
some structured sparsity in the data is present then one may even solve problems of 
potentially large size by using an appropriate sparse version of these LMI relaxations 
as described in [T7] , see also [HI §4.6]. 

The outline of the paper is as follows. In Section |2] we formally state the problem to be 
solved. In Section [3] we describe our hierarchy of inner approximations. In Section HJ we 
show that the specific geometry of the stability region can be exploited, as illustrated on 
several standard problems of robust control. The final section collects technical results 
and the proofs. 

2 Problem statement 

Let WL[x] denote the ring or real polynomials in the variables ), and let IRfxL 

be the vector space of real polynomials of degree at most d. Similarly, let E[x] C WL[x] 
denote the convex cone of real polynomials that are sums of squares (SOS) of polynomials, 
and £[x]d C S[x] its subcone of SOS polynomials of degree at most 2d. Denote by § m the 
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space ofmxm real symmetric matrices. For a given matrix A G § m , the notation A >z 
means that A is positive semidefinite, i.e., all its eigenvalues are real and nonnegative. 

Let P : IR[a;,tt] — > S m be a matrix polynomial, i.e. a matrix whose entries are scalar 
multivariate polynomials of the vector indeterminates x and u. Then 

P := {x G W 1 : Vw G U, P(x,u) h 0} (1) 

defines a parametrized polynomial matrix inequality (PMI) set, where x G MJ 1 is a vector 
of decision variables, u G M p is a vector of uncertain parameters belonging to a compact 
semialgebraic set 

U := {u G W : fli(u) > 0, z = 1, . . . , nj (2) 

described by given polynomials a,i(u) G M[-u], and P(x, u) is a given symmetric polynomial 
matrix of size m. As U is compact, without loss of generality we assume that for some 
i = i*, cii*(u) = R 2 — u T u, where R is sufficiently large. 

We also assume that P is bounded and that we are given a compact set B D P with 
explicitly known moments y = (y a ), ot G N n , of the Lebesgue measure on B, i.e. 

y a := / x a dx (3) 
Jb 

where x a := YYi=i X T ■ Typical choices for B are a box or a ball. To fix ideas, let 

B := {x G W 1 : bj{x) > 0, j = 1, . . . , n b } 

for some polynomials bj G Again, with no loss of generality, we may and will assume 
that for some j = j*, bj*(x) = R 2 — x T x, where R is sufficiently large. Finally, denote by 
vol A the Lebesgue volume of any Borel set A C B. 

We are now ready to state our polynomial inner approximation problem. 

Problem 1 (Inner Approximations) Given setP, build up a sequence of basic closed 
semialgebraic sets = {x G B : ga(x) > 0}, for some G M[x] } such that 

G d CP, d=l,2,... and lim vol G d = vol P. 

d— >oo 

In addition, we may want the sequence of inner approximations to satisfy additional 
nesting or convexity conditions. 

Problem 2 (Nested Inner Approximations) Solve Problem U\ with the additional 
constraint 

G d C G d+ i C P, rf=l,2,... 

Problem 3 (Convex Inner Approximations) Given set P, build up a sequence of 

nested basic closed convex semialgebraic sets Gd = {x G B : gd(x) > 0}, for some 
gd G M[x], such that 

G d C G d+l CP, d=l,2,... 
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3 A hierarchy of semialgebraic inner approximations 



Given a polynomial matrix P(x,u) which defines the set P in ([T]), polynomials a, G M[u] 
which define the uncertain set U in (J2J), let V = {v G lR m : v T v = 1} denote the Euclidean 
unit sphere of K m and let A m i n : B — > R be the function: 

x i — y A m ; n (a;) = min min v T P(x, u)v (4) 

as the robust minimum eigenvalue function of P(x,u). Function A m i n is continuous but 
not necessarily differentiable. It allows to define set P alternatively as the superlevel set 

P = {x G R™ : A min (s) > 0}. 



3.1 Primal SOS SDP problems 

Let ao G R[m] be the constant polynomial 1. Let 2d > max(2+degP, max, dega,, max,,- degbj), 
and consider the hierarchy of convex optimization problems indexed by the parameter 
d G N, d > d : 



Pd 



/ A min (x) dx — min / g{x) dx 

s.t. v T P(x, u)v — g(x) = r(x, u, v )(1 — v T v ) (5) 



+ Si(x, u, v)a,i{u) + tj(x, u, v)bj(x) V(x, u, v) 

i=0 j=l 

where decision variables are coefficients of polynomials g G M[x]2d, r G M>[x,u,v]2d r and 
coefficients of SOS polynomials Sj G T:[x,u,v]d s ., i = 0,1,..., n a , and tj G Ti[x,u,v]d t ., 
j = 1, ... ,71ft. Note in particular that the degrees of the polynomials should be such that 
d r > d — 1, d Si > d — |~(degai)/2] and d t > d — [(deg6j)/2"|. Since higher degree terms 
may cancel, the degrees can be chosen strictly greater than these lower bounds. However, 
in the experiments described later on in the paper, we systematically chose the lowest 
possible degrees. 

For each d 6 N fixed, the associated optimization problem (J5J) is a semidefinite pro- 
gramming (SDP) problem. Indeed, stating that the two polynomials in both sides of the 
equation in ([5]) are identical translates into linear equalities between the coefficients of 
polynomials g,r, (sj), (tj) and stating that some of them are SOS translates into semidef- 
initeness of appropriate symmetric matrices. For more details, the interested reader is 
referred to e.g. [HI Chapter 2]. 

3.2 Dual moment SDP problems 

To define the dual to SDP problem (J5]) we must introduce some notations. 
With a sequence y — (y a )> a G N n , let L y : M,[x] — > IR be the linear functional 

/ (=^/ a ^) h. L y (f) = J2f»y<*i / £R I4 
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With d G N, the moment matrix of order d associated with y is the real symmetric matrix 
M d (y) with rows and columns indexed in N d , and defined by 

M d (y)(a,(3) := L y (x a+ ?) = y a+p , Va,/3 G N n d . (6) 

A sequence y = (y a ) has a representing measure if there exists a finite Borel measure \i 
on M. n , such that y a = J x a dfi for every a G N n . 

With y as above and h G the localizing matrix of order d associated with y and h 
is the real symmetric matrix M d (hy) with rows and columns indexed by N d , and whose 
entry (a, (3) is given by 

M d {y){hy)(a,(3) := L^)^) = ^ /i 7 y a+/3+7 , Va, /3 G W d . (7) 
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With these notations, the dual to SDP problem (j3j) is given by: 

Pd = J Kam(x)dx — min L y (v T P(x, u)v ) 

s.t. M rf (y)hO, Afd_i((l-i; T T;)i/) = 



M d _ doi (a i2/ ) >r0, z = 0,l,...,n a ^ 
M d _ dbj (bjy) >z 0, j = l,...,n b 
Ly(x a ) = J B x a dx, Va G N" d 



where y G 



n+p+m 



3.3 Convergence 

Before stating our main results, let us recall some standard notions of functional analysis. 
Let g : B — > R be a function of a;, and let (<fa) denote a sequence of functions of x 
indexed by d G N. Lebesgue space ^i(B) is the Banach space of integrable functions on 
B equipped with the norm 




Regarding sequence (g d ), we use the following notions of convergence in B when d — > oo: 

• 9d — ► 9 in Li norm means lim \\g — g d \\i = 0; 

d— too 

• 9d g in Lebesgue measure means that for every e > 0, 

lim vol{x : — ^d(a?)| > s} = 0; 

• gd g almost everywhere means that lim^oo g d (x) = g(x) pointwise except possi- 
bly for x G A C B with vol A = 0; 

• 9d ~* 9 almost uniformly means that for every given e > 0, there is a set A C B 
such that vol A < e and gd^ g uniformly on B \ A; 
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• finally, with the notation g d f g we mean that g d — >■ (7 while satisfying ^(x) < 
5f d+ i(x) for all d. 

For more details on these related notions of convergence, see pfl §2.5]. 

Lemma 1 For every d > do, SDP problem has an optimal solution g d G R[x]2d and 

Pd = (>^mm(x) - g d (x)) dx = ||A min - g d \\i. (9) 
Jb 

A detailed proof of Lemma [I] can be found in §6.11 In particular we prove that there is no 
duality gap between SOS SDP problem ([5]) and moment SDP problem (jSJ), i.e. p d = p* d . 

For every d > do, let g d '■ B — > R be the piecewise polynomial 

:= max # fc (x). (10) 

do<k<d 

We are now in position to prove our main result. 

Theorem 1 Let g d G M[x]2d be an optimal solution of SDP problem and consider the 
associated sequence (gd) C -^i(B) for d > d . Then: 

(a) ^ — >■ A m ; n m Li norm and in Lebesgue measure; 

(b) g d t ^min almost everywhere, almost uniformly and in Lebesgue measure. 

A detailed proof of Theorem [I] can be found in §6.21 It relies on the Stone- Weierstrass 
theorem, Putinar's Positivstellensatz, Lebesgue's dominated convergence theorem and 
Egorov's theorem. 

3.4 Polynomial and piecewise polynomial inner approximations 

Corollary 1 For every d > do, let g d G IRfic^d be an optimal solution of SDP problem 
(EJ), let g d be the piecewise polynomial defined in ^^), and let 

G d := {x G B : g d {x) > 0}, G d := {x G B : g d {x) > 0}. (11) 

Then 

G d cP Vd>d and lim vol(P \ G d ) = 0. (12) 

G do C • • • C G d C • • • C P and lim vol(P \ G d ) = 0. (13) 

Tnat is, sequence (G d ) solves Problem^ and sequence {G d ) solves Problem® if piecewise 
polynomials are allowed. 

A proof can be found in §6.31 
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3.5 Nested polynomial inner approximations 



We now consider Problem [2] where gd is constrained to be a polynomial instead of a 
piecewise polynomial. We need to slightly modify SDP problem (jSJ). Suppose that at step 
d— 1 in the hierarchy we have already obtained an optimal solution g d _i G R[x] 2 d_2, such 
that gd-i > gd Q on B, for all do < d — 1. At step d we now solve SDP problem (J2]) with 
the additional constraint 

n b 

g(x) - gd-i{x) = cq(x) + ) j c j (x)b j (x), Vx (14) 

i=i 

with unknown SOS polynomials c G X[x] d and Cj G Efa^-^. . 

Corollary 2 Lei gr^ G R[x]2d be an optimal solution of SDP problem |3]) with the addi- 
tional constraint (fl^| ) and let Gd 6e as m < f77]j /or d > d . Then the sequence (G4) solves 
Problem^ 

For a proof see §6.41 

3.6 Convex nested polynomial inner approximations 

Finally, for g G Rfa^dj denote by V 2 g(x) the Hessian matrix of g at x, and consider SDP 
problem ([5} with the additional constraint 



for some SOS polynomials Co G £[x, v]d, Cj G E[x, i/^-^. anc ^ c «b+i e ^fo^ld-i- 

Corollary 3 Let g G R[x]2d be an optimal solution of SDP problem with the additional 
constraint (T73]) and let Gd be as in / TO]) for d > do- Then the sequence (Gd) solves Problem 



v T V 2 g{x)v = c Q (x,v) + ^2 / c j {x,v)b j (x) + c nb+1 (x,v)(l - v T v) 



(15) 



The proof follows along the same lines as the proof of Corollary [2j 



3.7 Example 



Consider the nonconvex planar PMI set 



P = {x G R 2 : P{x) 



1 — 16x1X2 



^0} 



,71 2 /y>2 

1 2 



which is Example II-E in [7j scaled to fit within the unit box 



B = {x G R 2 : 



A\oo < 1} 
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Figure 1: Degree two (left) and four (right) inner approximations (light gray) of PMI set 
(dark gray) embedded in unit box (dashed). 



whose moments ([3]) are readily given by 

4 



Va 



(ai + l)(a 2 + l)' 



On Figured] we represent the degree two and degree four solutions to SDP problem 
modelled by YALMIP 3 and solved by SeDuMi 1.3 under a Matlab environment. We see 
in particular that the degree four polynomial superlevel set G 2 is somewhat smaller than 
expected. This is due to the fact that the objective function in problem <^ is the integral 
of g(x) over the whole box B, not only over PMI set P. There is a significant role played 
by the components of the integral on complement set B\P, and this deteriorates the inner 
approximation. 

This issue can be addressed partly by embedding P in a tighter set B, for example here 
the unit disk 

B = {x e R 2 : \\x\\ 2 < 1} 
whose moments ([3]) are given by 



Vo 



r(so±i)r(2^i; 



r(2 



2 - 



where V is the gamma function such that T(k) = (k — 1)\ for integer k. See [121 Theorem 
3.1] for the general expression^] of moments of the unit disk in R™. 

On Figure [2] we represent the degree two and degree four solutions to SDP problem (jSJ). 
Comparing with Figure [TJ we see that the approximations embedded in the unit disk are 



2 Note however that there is an incorrect factor 2 ™ in the right handside of equation (3.3) in [12] . 
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Figure 2: Degree two (left) and four (right) inner approximations (light gray) of PMI set 
(dark gray) embedded in unit disk (dashed). 




x i 



Figure 3: Degree six (left) and eight (right) inner approximations (light gray) of PMI set 
(dark gray) embedded in unit disk (dashed). 

much tighter than the approximations embedded in the unit box. Finally, on Figure |3] 
we represent the tighter degree six and degree eight inner approximations within the unit 
disk. 
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4 Geometry of control problems 



As explained in the introduction, inner approximations of the stability regions are essen- 
tial for fixed-order controller design. The PMI regions arising from parametric stability 
conditions have a specific geometry that can be exploited to improve the convergence 
of the hierarchy of inner approximations. In this section, we first recall Hermite's PMI 
formulation of (discrete-time) stability conditions. Then we recall that the PMI stability 
region is the image of a unit box through a multi-afnne mapping, which allows to derive 
explicit expressions for the moments of the full-dimensional stability region, as well as 
tight polytopic outer approximations of low- dimensional affine sections of the stability re- 
gion. Numerical examples illustrate these techniques for fixed-order nominal and robustly 
stabilizing controller design. 



4.1 Hermite's PMI 



Derived by the French mathematician Charles Hermite in 1854, the Hermite matrix cri- 
terion is a symmetric version of the Routh-Hurwitz criterion for assessing stability of 
a polynomial. Originally it was derived for locating the roots of a polynomial in the 
open upper half of the complex plane, but with a fractional transform it can be readily 
transposed to the open unit disk and discrete-time stability. The criterion says that a 
polynomial x(z) = z n + x\z n ~ x + ■ • • + x n -\z + x n has all its roots in the open unit disk 
if and only if its Hermite matrix P{x) = Tj \x)Ti{x) — T^(x)T 2 (x) is positive definite, 
where 

1 Xi x 2 
1 xi 



T 2 (x) 



%n %n—l %n—2 
X n X n —\ 

x„ 



are n-by-n upper-right triangular Toeplitz matrices, see e.g. the entrywise formulas of 
[21 Theorem 3.13] or the construction explained in [I]. The Hermite matrix is n-by-n, 
symmetric and quadratic in coefficients x = (x\,X2, ■ ■ ■ ,x n ), so that the interior of the 
PMI set 

P = {i£l" : P(x) >z 0} 

is the parametric stability domain which is bounded, connected but nonconvex for n > 3. 
Optimal controller design amounts to optimizing over semialgebraic set P. 



4.2 Multiaffine mapping of the unit box 

As explained e.g. in [T5] or [TB^ §3.5] and references therein, stability domain P can also be 
constructed as the image of the unit box (in the space of so-called reflection coefficients) 
through a multiaffine mapping. More explicitly P = /(K) where K = {k 6 ~R n : ||fc||oo < 
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1} and multiaffine mapping / : M. n — > K n is denned by 



/(*) 









1 




h 







1 

h 
o 







h 
l 






h 
o 



1 





1 





A-2 





i + h 





h 





1 












l h 
fci l 
o o 



l 





k 2 k 3 + h(l + k 2 ) 
k 2 + kih(l + k 2 ) 
h 

in the case n = 3. The general expression of / for other values of n is not given here for 
space reasons, but it follows readily from the construction outlined above. 

Using this mapping we can obtain moments of B = P analytically: 

y a = f x a dx= [ (hh + hil + k^ih + hhil + k^k^detVfi^dk (16) 
Jp Jk 

where det V/(/c) = (1 + k 2 )(l — is the determinant of the Jacobian of /, in the case 
n = 3. For space reasons we do not give here the explicit value of function of a, 

but it can be obtained by integration by parts. 

Finally, let us mention a well-known geometric property of P: its convex hull is a polytope 
whose vertices correspond to the n + 1 polynomials with roots equal to —1 or +1. For 
example, when n = 3, we have 



convP = conv{(-3,3,-l), (-1,-1,1), (1,-1,-1), (3,3,1)}. 



(17) 



4.3 Third degree stability region 



Consider the problem of approximating from the inside the nonconvex stability region P 
of a discrete-time third degree polynomial z z 3 + X\z 2 + x 2 z + x 3 . An ellipsoidal inner 
approximation was proposed in |4J. The Hermite polynomial matrix defining P as in ([TJ 
is given by 



P(x) 



1-4 



Xi - X 2 X 3 



Xi 

x 2 



x 2 x 3 
Xlx 3 



l + x( 

X\ - 



x 2 

Xi 



x x x 3 

X 2 X 3 



X 2 X 3 



1-xl 



The boundary of P consists of two triangles and a hyperbolic paraboloid. The convex hull 
of P is the simplex described in ( TT7|) . We have analytic expressions (fl6|) for the moments 
© of B = P. 

On Figures IH [5] and |6] we respectively represent the degree two, four and six inner ap- 
proximations of P, scaled within the unit box for visualization purposes. We observe that 
the degree six approximation is very tight, thanks to the availability of the moments of 
the Lebesgue measure on P. 
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Figure 4: Two views of a degree two inner approximation (red) of nonconvex third-degree 
stability region (gray). 




Figure 5: Two views of a degree four inner approximation (red) of nonconvex third-degree 
stability region (gray). 
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Figure 6: Two views of a degree six inner approximation (red) of nonconvex third-degree 
stability region (gray). 



4.4 Fixed-order controller design 

Consider the linear discrete-time system with characteristic polynomial zhiz 4 - (2xi + 
X2)z i + 2x\z + 22 depending affinely on two real design parameters x\ and xi. It follows 
from Hermite's stability criterion that this polynomial has its roots in the open unit disk 
if and only if 

1 — x\ —2X\ — X2 — 2X\X2 2x\ + 2X\X2 + x\ 

„/ \ — 2x\ — X2 — 2x\X2 1 + 4xiX9 — 2x\ — xo — 2x\X2 
J (x) 

— 2x± — X2 — 2x\X2 1+4x1^2 — 2x\ — X2 — 2x\X2 

2X\ + 2X\X2 + x\ — 2X\ — X2 — 2X\X2 \ — x\ 

is positive definite. As recalled in (14. 2p . the convex hull of the four- dimensional sta- 
bility domain of a degree four polynomial is the simplex with vertices (—4,6,-4,1), 
(-2, 0, 2, -1), (0, -2, 0, 1), (2, 0, -2, -1), (4, 6, 4, 1) corresponding to the five polynomials 
with zeros equal to —1 or +1. Using elementary linear algebra, we find out that the 
image of this simplex through the affme mapping (— (2xi + X2), 0, 2xi, X2) parametrized 
by s 6 I 2 is the two-dimensional simplex 

B = co n v{(--, !),(-,--), (--,--)}. 

The (closure of the) stability region P = {x G M 2 : P(x) >z 0} is therefore included in B, 
whose moments ([3]) are readily obtained e.g. by the explicit formulas of [TO] . 
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Figure 7: Degree two (left) and four (right) inner approximations (light gray) of PMI 
stability region (dark gray) embedded in simplex (dashed). 
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Figure 8: Degree six (left) and eight (right) inner approximations (light gray) of PMI 
stability region (dark gray) embedded in simplex (dashed). 

On Figures and [8] we represent the degree two, four, six and eight inner approximations 
to P, corresponding to stability regions for the linear system. We observe that the ap- 
proximations become tight rather quickly. This is due to the fact that B is a good outer 
approximation of P with known moments. Tighter outer approximations B would result 
in tighter inner approximations of P, but then the moments of B can be hard to compute, 
see |8]. 



4.5 Robust controller design 

Now consider the uncertain polynomial z \- > x<i + u + 2x\z — (2x\ + X2)z 3 + z 4 with 
!i 6 U = {ti 6 1 : M 2 < |} with uncertain Hermite matrix P(x, u) and the corresponding 
parametrized PMI stability region P in (pQ). Let us use the same bounding set B as in 

S3 
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Figure 9: Degree two (left) and four (right) inner approximations (light gray) of robust 
PMI stability region (dark gray) embedded in simplex (dashed). 

On Figure [9] we represent the degree two and degree four inner approximations to P, 
corresponding to robust stability regions for the linear system. Comparing with Figure 
[7] we see that the approximations are smaller, and in particular they do not touch the 
stability boundary to cope with the robustness requirements. 

5 Conclusion 

We have constructed a hierarchy of inner approximations of feasible sets defined by 
parametrized or uncertain polynomial matrix inequalities (PMI). Each inner approxi- 
mation is computed by solving a convex linear matrix inequality (LMI) problem. The 
hierarchy converges in a well-defined analytic sense, so that conservatism of the approx- 
imation is guaranteed to vanish asymptotically. In addition, the inner approximations 
are simple polynomial or piecewise-polynomial superlevel sets, so that optimization over 
these sets is significantly simpler than optimization over the original parametrized PMI 
set. In particular, we remove the possibly complicated dependence of the problem data 
on the uncertain parameters. 

One may also impose the hierarchy of inner approximations to be nested. Finally, one may 
also impose the inner approximations to be convex. In this latter case they do not converge 
any more to the feasible set but, on the other hand, optimization over the parametrized 
PMI set can be reformulated as a convex polynomial optimization problem (of course at 
the price of some conservatism). Ideally, beyond convexity, we may also want the inner 
convex approximation to be semidefinite representable (as an explicit affine projection 
of an affine section of the SDP cone), and deriving such a representation may be an 
interesting research direction. 

The tradeoff to be found is between tightness of the inner approximation and degree of 
the defining polynomials. In the context of robust control design, a satisfactory inner 
approximation can be possibly computed off-line, and then used afterwards on-line in a 
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feedback control setup. 

Our methodology is valid for general parametrized PMI problems. However, in the case of 
parametrized PMI problems coming from fixed-order robust controller design problems, 
geometric insight can be exploited to improve convergence of the hierarchy. The key 
information is the knowledge of the moments of the Lebesgue measure on a compact 
set which tightly contains the parametrized PMI set we want to approximate from the 
inside. In turns out that for robust control problems this knowledge is available easily, as 
illustrated in the paper by several examples. 

The main limitation of the approach lies in the ability of solving primal moment and dual 
polynomial sum-of-squares LMI problems. State-of-the-art general-purpose semidefinite 
programming solvers can currently address problems of relatively moderate dimensions, 
but problem structure and data sparsity can be exploited for larger problems. 

Acknowledgements 

The first author acknowledges support by project No. 103/10/0628 of the Grant Agency 
of the Czech Republic. We are grateful to Luca Zaccarian for pointing out a mistake in a 
previous version of this paper. 

6 Appendix 

6.1 Proof of Lemma [1] 

Proof: The dual of polynomial SOS SDP problem (j5J) is moment SDP problem (181 . 
Slater's condition cannot hold for (jHj) because V has empty interior in M. m . However it 
turns out that Slater's condition holds for an equivalent version of SDP problem fl8]), i.e., 
the latter has a strictly feasible solution y. Indeed, let J C M[v] be the ideal generated by 
the polynomial v i— > 9{v) := 1 — v T v so that the real variety Vk(J) := {v G M m : 6(v) = 0} 
associated with J is just the unit sphere V. It turns out that the real radical of J is J 
itself, that is, /(Vr(J)) = J (where for S C M m , I(S) denotes the vanishing ideal). And 
after embedding J in u, v], we still have I(Vr(J)) = J. 

Let H := {(ot,/3, 7 ) G N n x W x N m : 7m < 1}, and let H d := {(a,/3, 7 ) G H : 
J2i a i + Yuj Pj + 7^ — T ne monomials (x a u^v 7 ), (a, (3, 7) G H, form a basis of the 
quotient space M[x,u,v]/J. Moreover, for every (at, /3, 7) G pj" +p+m ; 

X° L V?V 1 = 2j PabcX a U b V C + h(x,U,v) (1 - V T v), 
(aAc)eHd &[x,u,v\d-2 

for some real coefficients (jw), and some h G M.[x, u, u]d_2- This is because every time one 
sees a monomial x a u b v c with c m > 2, one uses v ^ = 1 — X^y m v j ^° re duce this monomial 

3 V is Zariski dense in V C (J) (= {v G C™ : 0(v) = 0}) so that I(V R {J)) = I{V C (J)). But 6 being 
irreducible, J is a prime ideal and so I(Vc(J)) = J. 
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modulo 9 = (1 — v T v). For instance 

x a u b v?---v c - 1 1 v 3 m = AV-C"i\x v 2 m 



x°uV Cl • • • o^d™ - V xV*;? 1 • • • v? +2 ■ ■ ■ v c Sl'v 



1 ^ m _i ^ m ^ ^ u, u x uj "m-1 v ™ 

V v ' 



eM[x,u,v]d-2 

etc. Therefore, for every (a,/3, 7), (a', /?', 7') G i^d, 



2j p abc x a u b v c + h(x, u,v)^ (1 - f T w), 

(a,&,c)€#2d 



for some real coefficients (p a bc), and some h G Rfsc, w, f]2d-2- 

So because of the constraints M ( j_i((l — v T v)y) = 0, the semidefinite program OH]) is 
equivalent to the semidefinite program: 



Pi 



*d= 1 ^mm(x)dx — min L y (v T P(x,u)v) 

Jb y 

s.t. M d {y) h 0, M d _!((l - v T i;) y) = 
Mi-d aj (a. y) y 0, z = 0,l,...,n o 
M d -d bj {bjV) b.0, j = l,...,n b 
L y {x a ) = J B x a dx, VaeN^ 



(19) 



where the smaller moment matrix M d (y) is the submatrix of M d {y) obtained by looking 
only at rows and columns indexed in the monomial basis (x a y l3 v~ t ), (a,/3,j) G H d , in- 
stead of f^™ +p+m . Similarly, the smaller localizing matrix M d ^ da .(aiy) * s the submatrix 
of Md-da-i^iy) obtained by looking only at rows and columns indexed in the monomial 
basis (x^y^f 7 ), (a, /3, 7) G H d _ da ., instead of N^^ +m ; and similarly for M d _ d (b-y). 

Indeed, in view of ( TT8|) and using M d ((l — v T v) y) = 0, every column of M d (y) associated 
with (a,/3, 7) G pj ra +P +?n is a linear combination of columns associated with (a',/3',7') G 
H d . And similary for M d _ da ia iy ) and M d -d b .(b jy ). Hence, M d (y) b <^ M d (y) b 0, 
and 

M d _ da . (a, y) y <S> M d _ d „. (a< y) b 0; M d _ 4 . (6 i y) b ^ M^.tyj !/)b0, 

1 1 J J 

for alH = 1, . . . , n a , j = 1, . . . , n&. 

Next, let y be the sequence of moments of the (product) measure \i uniformly distributed 
on B x U x V, and scaled so that for all (a, f3, 7) G N^J"^" 1 

y a/ 3 7 = A x a u^v' 7 dfi(x,u,v) = — — — — III x a v 1 dx du d\{v) 
JBxUxV volU x V J B Ju J v v v ' 

d/j,(x,u,v) 
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(with A the rotation invariant measure on V). Therefore, for every a G N^, 



y a QQ = L y (x a ) = / x a dp(x, u, v ) = / x a dx. 

JbxUxV Jb 

Moreover, M d _ 1 ((l—v T v) y) = for every d and importantly, M d (y) >~ 0, M d -d a . ( a « 2/) ^ 
and Md-d b . ipj y) >~ 0. To see why, suppose for instance that h T M d (y)h = for some vector 
h 7^ 0. This means that for some non trivial polynomial h G Wi[x,u,v]/ J of degree d, 



h T M d (y)h = h 2 dp = 0, 

JbxUxV 



that is, h(x, u,v) =0 for //-almost all (x, u, v) G B x U x V, and so h(x, u, v ) = for 
all (x, -u, f ) G B x U x V because is continuous. But as B x U has nonempty interior 
in M n x R p , then necessarily h G J(Vr(J)) (= J) - see Lemma [2] in section RTol - which 
contradicts ^ h G K[x, it, u]/ J. Therefore y is a strictly feasible solution of (JTHJ) and so 
Slater's condition holds for ffT9l). 



Denote by S[x, v] d the space of polynomials of degree at most 2c?, that are SOS of 
polynomials in M[x, u, v]/J. As y is a strictly feasible solution of the semidefinite program 
(j!9p . by a standard result of convex optimization, there is no duality gap between (T3~9|) 
and its dual 



p' d = A m i n (x) da; — min / g(a;) dx 



s.t. v T P(x, u)v — g(x) = r(x, u, v )(1 — v T v ) 

+ Si (a;, m, v)a,i(u) + £j(£> w> v)bj(x) V(x, u, u) 

i=0 j=l 

(20) 

where now the decision variables are coefficients of polynomials g G Kjx^d, r G w, f ]2d r , 
and coefficients of SOS polynomials s« G £[x, w, f]d a -> « = 0, 1, . . . , n a , and G u, v] db ,, 
j = 1, .. .,7ii,. That is, = and so = p* d because p' d < p d < p* d . If p' d < oo then 
(120]) is guaranteed to have an optimal solution (g*,r*,s*,t*). But observe that such an 
optimal solution (g*,r*,s*,t*) is also feasible in ([5]), and so having value p' d = p d = p d , 
(g*, r*,s*, t*) is also an optimal solution of (JSJ). 

It remains to prove that p d is bounded. For any feasible solution y of (J5]), yo < volB, and 

L y (x?) < R 2d y d o; L y {uf) < R 2d y$\ Lyttf) < vt (21) 

for all i = 1, . . . , n, j = 1, . . . ,p, k = 1, . . . ,m. This follows from M d _ da t (a«. y) >z 0, 
M d _ db t (bj* y) >z and Md_i((l — v T v)y) = 0, where cii*(x) = R 2 — x T x and bj*(x) = 

R 2 — u T u; see the comments after (T5]) and ([3]). Then by [131 Lemma 4.3], one obtains 
\Ua\ < i? 2d (volB) d , for all a G N^, which shows that the feasible set of (JHJ) is compact. 
Hence (jSJ) has an optimal solution and p d is finite; therefore its dual flSJ) also has an optimal 
solution, the desired result. □ 
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6.2 Proof of Theorem [T] 

Proof: (a) Let K:=BxUxVc M n +P +m and consider the infinite-dimensional opti- 
mization problem 

p = min / v T P(x,u)v dp(x,u,v) 

M eM(K) J K ^ 

s.t. / x a dp = [ x a dx, a G N n 



K 



where M(K) is the space of finite Borel measures on K. Problem ( 122]) has an opti- 
mal solution /i* G M(K). Indeed, p > f B \ min (x)dx because for every (x, u,v) G K, 
v T P(x,u)v > A min (:r); and so for every feasible solution p G M(K), 



v P(x,u,v)v dp(x,u,v) > / \ m i n (x) dfi(x,u,v) = / A m i n (x) ofe 
k Jk Jb 

because J K x a dp = J B x a dx for all a G N and hence the marginal of p on W 1 is the 
Lebesgue measure on B. On the other hand, observe that for every x G B, A m i n (x) = 
v^P(x, u x )v x for some (u x , v x ) G U x V. Therefore, let p* G M(K) be the Borel measure 
concentrated on (x,u x ,v x ) for all x G B, i.e. 

//(B' x U' x V) := f lu/ xV / ( Ua ., v.,) dx, V(B', U', V) G B(B) x B(U) x B(V) 

JB'nv 

where x h- >■ 1b(x) denotes the indicator function of set B and -B(B) denotes the Borel 
a- algebra of subsets of B. Then p* is feasible for problem (|22p with value 



v P(x,u)v dp*(x,u,v) = / X m i n (x) dx 
k Jb 

which proves that p = f B A m i n (x) dx. 

Next, A min being continuous on compact set B, by the Stone-Weierstrass theorem [TJ 
§A7.5], for every e > there exists a polynomial h £ G M[x] such that 

sup |A min (x) - h e (x)\ < |. 

x&B * 

Hence the polynomial p £ := h e — e satisfies A min — p e > on B and so v T P(x, u)v — p £ > 
on B x U x V. By Putinar's Positivstellensatz, see e.g (TTj Section 2.5], there exists SOS 
polynomials r e G M[x,M,f], and Si e ,tj £ G T,[x,u,v] such that equation (j^J) is satisfied. 
Hence for d sufficiently large, say d > d £ , (p £ ,r e , Si e ,tj £ ) is a feasible solution of (E]) with 
associated value 

J (Kmn(x) -p £ {x))dx < y J dx. 

Hence < pa < y J s ofe whenever d > d £ where is defined in (j§J) • As s > was 
arbitrary, we obtain the desired result 

lim p d = 0. 

d— too 
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Observe that since g d < A min for all d, 

Pd= (Amin (a;) - 9d(x)) dx = I |A min (a;) - g d (x)\ dx 



B 



so that the convergence p d — > is just the convergence g d —> A min for the L\ norm on B. 
Finally the convergence g d — > A mm in Lebesgue measure on B follows from [TJ Theorem 
2.5.1]. 

(b) For each x G B, fixed and arbitrary, the sequence (g d ) is monotone nondecreasing and 
bounded above by A mm . Therefore there exists g* : B — > R such that for every x G B, 
^d(^) t 5 l *( a; ) < A mm (a;) as <i — > oo. Since ^ > #o an d j B godx > — oo, by Lebesgue's 
Dominated Convergence Theorem P, §1.6.9] 

g*(x)dx = lim / g d (x)dx = / A mm (a;)(ia;, 

B Jb Jb 

and so from #*(x) < A min (:r) we deduce that g*(x) = A min (x) for almost all x G B. 
Combining the latter with g d \ g*, we obtain that g d — > A m i n almost everywhere in B. 
But then since the Lebesgue measure is finite on B, by Egorov's theorem (TJ Theorem 
2.5.5], g d A m i n almost uniformly in B. Finally, convergence in Lebesgue measure on B 
also follows from [TJ Theorem 2.5.2]. □ 



6.3 Proof of Corollary [T] 

Proof: By Theorem [TJ lim^oo ||A min — g d \\i = 0. Therefore, by [TJ, Theorem 2.5.1] the 
sequence (go) converges to A min in Lebesgue measure, i.e. for every e > 0, 

lim vol{x : |A min (x) - g d {x)\ > e} = 0. (23) 

d— >oo 

Let £ > be fixed, arbitrary, and let P e := {x G B : A mm (x) > e}, so that lim £ ^ vol P £ = 
volP. By §25$), lim d _, 00 vol(P e n {x G B : g d (x) < 0}) = 0. Next, for all d G N, 

volP £ = vol(P £ n {x G B : g d {x) < 0}) + vol(P £ n {x G B : g d (x) > 0}). 

Therefore, taking the limit as d — > oo yields 

volP e = lim vol(P e n {x G B : g d (x) < 0}) + lim vol(P £ n {x G B : g d (x) > 0}) 

=o by (H3J) =Gd 
= lim vol(P e n G d ) < lim vol G d . 

As e > was arbitrary and C P, we obtain the desired result ( TT21 . The proof of (TT5|) 
is similar. □ 
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6.4 Proof of Corollary H 



Proof: Let < e < | be fixed, arbitrary. As in the proof of Theorem [H for every fcGN 
there exists a polynomial h k G M[x] such that sup.j, GB |A min (:c) — hk(x)\ < e k . Hence for 
all x G B and all k > 1, 

X min (x)-3e k < h k (x)-2e k < \ min (x)-s k < \ min (x)-3s k+1 < h k+1 (x)~2e k+1 < A min (x)-e 

and so the polynomial x p k (x) := h k (x) — 2e k satisfies p k+ i(x) > p k (x) and A m i n (x) > 
p k (x) for all x G B. Again, by Putinar's Positivstellensatz, see e.g [H] Section 2.5], p k is 
feasible for with the additional constraint (THj) . provided that is sufficiently large, 
and with associated value 



k+l 



|A m i n (x) — p k (x)\dx = I (X mm (x) — p k (x))dx < 3e \ dx — > as k — > oo. 

□ 



6.5 An auxiliary result for the proof of Lemma [T] 

Remember that J C M[v] is the ideal generated by 1 — v T v and the real radical I(Vr(J)) 
of J is J itself. And when J is embedded in Wi[x, u, v] (with same name of simplicity) we 
also have /(Vr(J)) = J. 

Lemma 2 If f £ M[x, u, v]d is such that f(x, u, v) — for all (x, u, v) G B x U x V then 

feJ. 

Proof: Write 

f(x,U,v) = ^2 9a(x,u)v a , 

for some polynomials (g Q ) C tt]d, a G N™. Next, let (x ,u ) G B x U be fixed, so 
that v i — y f(xo, uq, v ) = for all v G V. Therefore, as a polynomial of R[i>], it vanishes on 
V = Vr(J) and as I(Vr(J)) = J, v /(xo,Mo,w) G J, that is, 

/(z ,u ,t;) = 9a(x ,u )v a = (l-v T v)9 x °' v °{v), (24) 

for some polynomial t> i— >■ fl 20 '" ^) G K[f]d- The coefficients (#^ ' u °) of the polynomial 
Qx ,u ^ _ ^ Q ^0'"0ti a are linear in the coefficients (gp(x , w )), /3 G N™, of /. Indeed 
one may reduce each monomial v a using t>^ = 1 — ^^m^ 2 ' un t n there is no monomial 
with /3 > 1. For instance, 

vT ■ ■ ■ vT-xvl = v? ■ ■ ■ v^ 1 (v T v - 1) + v? ■ ■ ■ - £ v\ 
and 



otj+2 



i-l 

-1 5 
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etc., to finally obtain 

v a =p a (v)(l-v T v)+r a , Va G N™, 
for some p a G R[w]i a |_ 2 and r a G M[v]/J. Therefore, summing up over all a G N™ yields: 

f(x Q ,uo,v) = ^2 g a (x ,u )v a 

= (v T v-l) ^ g a {x ,uo)p a (v) + ^ g a (xo,u ) r a (v) 



*GN™ ct€N™ 



ft,(x ,«o,t)) = o as u h-> /(x , Mo, f ) G J 

= (w T w- l)/i(a;o,uo,w), VwGR m , (25) 

for some h G R[x,M,f]. But since (125]) holds for every (x,u) G B x U, we obtain 

f(x,u,v) = (v T v - 1) h(x ,Uo,v), V(x, u, v ) G B x U x ]R m , 

and as B x U x R m has nonempty interior, 

f(x, u, v) = {v T v - 1) h(x , M , u), V(x, w, u) G M n x R p x M m , 

i.e., / = {v T v — l)h, which proves the desired result that / G J. □ 
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